/// 基数は10の冪の多倍長整数。
/// 基数以下の配列で絶対値（負でない数）を表す。0は空配列。先頭にムダな0は許さない。
/// 一度作った絶対値は原則破壊されない（複数のBignumから参照可能。演算は遅いがコピーは一瞬）
/// 絶対値はリリース版ではチェックされない（明示的にabschkしない限り）。直接作らない限り無問題のはず。
/// Date: September 8, 2008
module bignum;
import std.string;

private int pow10(int n) {
  int ret=1;
  for( ; n>0; n--) ret*=10;
  return ret;
}

/// N=10^Kが基数となる。1≤K≤9で、2N+1がintの範囲に収まる。
const K=9;
static assert(1<=K && K<=9);
const N=pow10(K); /// ditto

/// ゼロで除算したとき投げる
class ZeroDivException : Exception {
  this() { super("Division by zero!"); }
}

/// 絶対値への変換失敗。
class AbsConvError : Error {
  this() { super("cannot convert to Absolute!"); }
}

/// Bignumに変換に失敗
class BigConvError : Error {
  this() { super("cannot convert to Bignum!"); }
}

/// 整数型への変換に失敗
class IntConvError : Error {
  this() { super("cannot convert to Integer!"); }
}

/// Bignumを整数型Tに変換
/// Bugs: Tは(u)intと(u)longのみ
/// Throw: Tの範囲外のときIntConvError
T toint(T)(Bignum b) {
  static Bignum max, min;
  if(!max) max=tobig(T.max);
  if(!min) min=tobig(T.min);

  if(bigcomp(b,min)<0 || bigcomp(max,b)<0) {
    throw new IntConvError();
  }

  T ret;
  for(size_t i=b.abs.length; i>0; i--) {
    ret = ret*N + b.sign*b.abs[i-1];
  }
  return ret;
}

/// new Bignum(s)と同じ
/// Throw: 変換失敗時BigConvError
Bignum tobig(string s) { return new Bignum(s); }

/// new Bignum(n)と同じ
Bignum tobig(int n) { return new Bignum(n); }
Bignum tobig(uint n) { return new Bignum(n); } ///ditto
Bignum tobig(long n) { return new Bignum(n); } ///ditto
Bignum tobig(ulong n) { return new Bignum(n); } ///ditto
Bignum tobig(int[] n) { return new Bignum(n); } ///ditto

/// そのまま返す
Bignum tobig(Bignum b) { return b; }

/// 数字列を絶対値に変換。先頭のゼロは無視。
/// Throw: 変換失敗時、AbsConvError
int[] toabs(string s) {
  if(s.length==0) throw new AbsConvError();

  int i;
  while(i<s.length && s[i]=='0') i++;

  if(s.length-i!=countchars(s[i..$],"0-9")) throw new AbsConvError();

  int[] a;
  int j=s.length;

  while(j-K>=i) {
    a~=atoi(s[j-K..j]);
    j-=K;
  }
  if(i<j) a~=atoi(s[i..j]);

  assert(i<s.length || a.length==0);
  return a;
}

/// 整数型から絶対値に変換
/// Throw: 負の数与えたときAbsConvError
int[] toabs(uint n) {
  if(n==0) return [];

  int[] abs;
  while(n>=N) {
    abs~=cast(int)(n%N);
    n/=N;
  }
  abs~=cast(int)n;

  return abs;
}

/// ditto
int[] toabs(ulong n) {
  if(n==0) return [];

  int[] abs;
  while(n>=N) {
    abs~=cast(int)(n%N);
    n/=N;
  }
  abs~=cast(int)n;

  return abs;
}

/// ditto
int[] toabs(int n) {
  if(n<0) throw new AbsConvError();
  return toabs(cast(uint)n);
}

/// ditto
int[] toabs(long n) {
  if(n<0) throw new AbsConvError();
  return toabs(cast(ulong)n);
}

/// 絶対値の正当性チェック。assert(abschk(abs))のように使う
bool abschk(int[] abs) {
  foreach(i; abs) {
    if(!(i>=0 && i<N)) return false;
  }
  return abs.length==0 || abs[$-1]>0;
}

/// 絶対値を文字列に変換。
string tos(int[] abs) {
  assert(abschk(abs));

  if(abs.length==0) return "0";

  char[] s;
  static string fom;
  if(!fom) fom=format("%%0%dd",K);

  for(size_t i=0; i<abs.length-1; i++) {
    s=format(fom,abs[i])~s;
  }
  s=format("%d",abs[$-1])~s;

  return cast(string)s;
}

/// 破壊的インクリメント(カウンタにどうぞ)
int[] absinc(int[] n) {
  size_t i=0;
  while(i<n.length && n[i]==N-1) n[i]=0;
  if(i<n.length) {
    n[i]+=1;
  } else {
    n~=1;
  }
  return n;
}

// 先頭ゼロ取り除く
private int[] delzero(int[] n) {
  size_t i=n.length;
  while(i>0 && n[i-1]==0) i--;
  return n[0..i];
}

///多倍長整数クラス
class Bignum {
  int sign;  /// 符号。0か±1。
  int[] abs; /// 絶対値。書き換えないこと（中身も）⇒D2.0ならinvariant?

  invariant() {
    assert(abschk(abs));
    if(abs.length>0) {
      assert(sign==1 || sign==-1);
    } else {
      assert(sign==0);
    }
  }

  private this() {}

  /// signとabsを直接初期化（ノーチェック）
  this(int s, int[] a) {
    sign=s;
    abs=a;
  }

  /// int[]で（非負に）初期化（チェックなし）
  this(int[] a) {
    this(abs.length ? 1 : 0, a);
  }

  /// 整数で初期化。
  this(uint n) {
    sign = n>0 ? 1 : 0;
    abs = toabs(n);
  }

  /// ditto
  this(ulong n) {
    sign = n>0 ? 1 : 0;
    abs = toabs(n);
  }

  /// ditto
  this(int n) {
    if(n>=0) {
      this(cast(uint)n);
    } else {
      this();
      sign=-1;
      abs=toabs(-cast(long)n);
    }
  }

  /// ditto
  this(long n) {
    static assert(-long.max-1==long.min);

    if(n>=0) {
      this(cast(ulong)n);
    } else {
      this();
      sign=-1;
      abs = n!=long.min ? toabs(cast(ulong)-n) : toabs(cast(ulong)long.max+1);
    }
  }

  /// 文字列で初期化。符号があってもよい（１つだけ）。符号以下はtoabsを使う。
  /// Throw: 変換失敗時BigConvError
  this(string s) {
    if(s.length==0) throw new BigConvError();

    try {
      switch(s[0]) {
        case '-':
          sign=-1;
          abs=toabs(s[1..$]);
          break;
        case '+':
          sign=1;
          abs=toabs(s[1..$]);
          break;
        default:
          sign=1;
          abs=toabs(s);
      }
    } catch(AbsConvError e) {
      throw new BigConvError();
    }

    if(abs.length==0) sign=0;
  }

  /// 文字列化メソッド
  override string toString() {
    return sign==-1 ? "-"~tos(abs) : tos(abs);
  }

  hash_t toHash() {
    return abs.length>0 ? abs[0] : 0;
  }

  int opEquals(Object o) {
    Bignum b = cast(Bignum)o;
    return bigcomp(this, b)==0;
  }

  int opCmp(Object o) {
    Bignum b = cast(Bignum)o;
    return bigcomp(this, b);
  }
}

/// 大小比較。n>mなら1、n==mなら0、n<mなら-1を返す。
int bigcomp(Bignum n, Bignum m) {
  if(n.sign>m.sign) return 1;
  if(n.sign<m.sign) return -1;

  return n.sign>=0 ? abscomp(n.abs,m.abs) : abscomp(m.abs,n.abs);
}

///ditto
int abscomp(int[] n, int[] m) {
  assert(abschk(n));
  assert(abschk(m));

  size_t ln=n.length, lm=m.length;

  if(ln>lm) return 1;
  if(lm>ln) return -1;

  while(ln>0 && n[ln-1]==m[ln-1]) ln--;
  if(ln==0) return 0;

  if(n[ln-1]>m[ln-1]) {
    return 1;
  } else {
    assert(m[ln-1]>n[ln-1]);
    return -1;
  }
}

/// 符号反転
Bignum bigminus(Bignum n) {
  return new Bignum(-n.sign,n.abs);
}

/// 足し算
Bignum bigadd(Bignum n,Bignum m) {
  if(n.sign==m.sign) {
    return new Bignum(n.sign,absadd(n.abs,m.abs));
  } else if(n.sign>m.sign) {
    return abssub(n.abs,m.abs);
  } else {
    return abssub(m.abs,n.abs);
  }
}

/// ditto
int[] absadd(int[] n, int[] m) {
  assert(abschk(n));
  assert(abschk(m));

  // nの方がケタ多い
  int[] absadd2(int[] n, int[] m) {
    assert(m.length<=n.length);

    if(m.length==0) return n;

    int[] c;
    size_t ln=n.length, lm=m.length;
    c.length=ln;

    int r=0;
    size_t i=0;
    for( ; i<lm; i++) {
      c[i]=n[i]+m[i]+r;
      r=(c[i]>=N);
      if(r) c[i]-=N;
    }

    for( ; r && i<ln; i++) {
      c[i]=n[i]+1;
      r=(c[i]==N);
      if(r) c[i]=0;
    }

    if(r) {
      assert(i==ln);
      c~=1;
    } else {
      c[i..$]=n[i..$];
    }

    return c;
  }

  if(n.length>=m.length) {
    return absadd2(n,m);
  } else { 
    return absadd2(m,n);
  }
}

/// 引き算。絶対値同士でも負になるかもしれないのでBignum返す。
Bignum bigsub(Bignum n, Bignum m) {
  return bigadd(n, bigminus(m));
}

/// ditto
Bignum abssub(int[] n, int[] m) {
  assert(abschk(n));
  assert(abschk(m));

  // nの方が大きい
  int[] abssub2(int[] n, int[] m) {
    assert(abschk(n));
    assert(abschk(m));
    assert(abscomp(n,m)>=0);

    size_t ln=n.length, lm=m.length;
    if(lm==0) return n;

    int[] c;
    c.length=ln;

    int r=0, max=-1;
    size_t i=0;

    for( ; i<lm; i++) {
      c[i]=n[i]-m[i]-r;
      r=(c[i]<0);
      if(r) c[i]+=N;
      assert(c[i]>=0);
      if(c[i]>0) max=i;
    }

    for( ; r; i++) {
      r=(n[i]==0);
      c[i] = r ? N-1 : n[i]-1;
      if(c[i]>0) max=i;
    }

    if(i<ln) {
      c[i..ln]=n[i..ln];
    } else {
      c.length=max+1;
    }

    return c;
  }

  size_t ln=n.length, lm=m.length;
  if(ln>lm) {
    return new Bignum(1,abssub2(n,m));
  } else if(lm>ln) {
    return new Bignum(-1,abssub2(m,n));
  }

  while(ln>0 && n[ln-1]==m[ln-1]) ln--;
  if(ln==0) return tobig(0);

  if(n[ln-1]>m[ln-1]) {
    return new Bignum(1,abssub2(n[0..ln], delzero(m[0..ln])));
  } else {
    return new Bignum(-1,abssub2(m[0..ln], delzero(n[0..ln])));
  }
}

/// 掛け算
Bignum bigmul(Bignum n, Bignum m) {
  int[] abs=absmul(n.abs, m.abs);

  if(abs.length>0) {
    int sign = n.sign==m.sign ? 1 : -1;
    return new Bignum(sign,abs);
  } else {
    return tobig(0);
  }
}

/// ditto
int[] absmul(int[] n, int[] m) {
  assert(abschk(n));
  assert(abschk(m));

  bool all0(int[] a) {
    foreach(i ; a) {
      if(i!=0) return false;
    }
    return true;
  }

  size_t ln=n.length, lm=m.length;
  if(ln==0 || lm==0) return [];

  int[] c;
  c.length=ln+lm;

  for(int i=0; i<lm; i++) {
    int r=0;

    for(int j=0; j<ln; j++) {
      long t = cast(long)n[j] * m[i] + r;
      assert(t <= cast(long)N*N-N);

      c[i+j]+=cast(int)(t%N);

      if(c[i+j]>=N) {
        c[i+j]-=N;
        r=cast(int)(t/N)+1;
      } else {
        r=cast(int)(t/N);
      }
      assert(r<N);
    }

    assert(all0(c[ln+i..$]));
    c[ln+i]=r;
  }

  if(c[$-1]==0) c.length=c.length-1;
  assert(c[$-1]>0);

  return c;
}

/// 割り算と余り
/// Throw: mが0のときZeroDivException
int[] absdiv(int[] n, int[] m) {
  int[] r;
  return absdivmod(n,m,r);
}

/// ditto
int[] absmod(int[] n, int[] m) {
  int[] r;
  absdivmod(n,m,r);
  return r;
}

/// ditto
Bignum bigdiv(Bignum n, Bignum m) {
  int[] abs=absdiv(n.abs, m.abs);

  if(abs.length>0) {
    int sign = n.sign==m.sign ? 1 : -1;
    return new Bignum(sign,abs);
  } else {
    return tobig(0);
  }
}

/// ditto
Bignum bigmod(Bignum n, Bignum m) {
  int[] abs=absmod(n.abs, m.abs);

  if(abs.length>0) {
    int sign = n.sign==m.sign ? 1 : -1;
    return new Bignum(sign,abs);
  } else {
    return tobig(0);
  }
}

/// 商を返すと同時にremに余りを入れる。
/// Throw: mが0のときZeroDivException
int[] absdivmod(int[] n, int[] m, out int[] rem) {
  assert(abschk(n));
  assert(abschk(m));

  //商<Nの割算 (m>0), rは余り
  int divN(int[] n, int[] m, out int[] r) {
    assert(abschk(n));
    assert(abschk(m));
    assert(m.length>0);
    assert(abscomp(n,[0]~m)<0);

    //3ケタ÷2ケタ(≤N)
    int div32(int[] a, int[] b) {
      assert(a.length==3 && b.length==2);
      assert(a[2]>0 && b[1]>0 && a[1]>=0 && a[0]>=0 && b[0]>=0);
      assert(a[2]<N && a[1]<N && a[0]<N && b[1]<N && b[0]<N);
      assert(a[2]<b[1] || (a[2]==b[1] && a[1]<=b[0]));

      long la=(cast(long)a[2]*N+a[1])/(b[1]+1);
      assert(0<la && la<N);
      int A=cast(int)la;

      Bignum bb = abssub(a,absmul([A],b));
      assert(bb.sign==1 && abscomp(bb.abs,[0,N-1,1])<0);

      long B;
      for(int i=bb.abs.length-1; i>=0; i--) B=B*N+bb.abs[i];
      long lb=cast(long)b[1]*N+b[0];

      int ret=A + cast(int)(B/lb);
      assert(ret<N || (a[2]==b[1] && a[1]==b[0] && ret==N));
      return ret;
    }

    size_t ln=n.length, lm=m.length;
    if(lm>ln) { r=n; return 0; }

    assert(lm==ln || lm+1==ln);
    if(ln==1) {
      r = n[0]%m[0] ? [n[0]%m[0]] : [];
      return n[0]/m[0];
    }

    if(ln==2) {
      long l = cast(long)n[1]*N+n[0];
      long k = lm==2 ? cast(long)m[1]*N+m[0] : m[0];
      long lr = l%k;
      if(lr>=N) {
        r = [cast(int)(lr%N),cast(int)(lr/N)];
      } else {
        r = lr ? [cast(int)lr] : [];
      }
      return cast(int)(l/k);
    }

    int p;
    if(lm==ln) {
      long k=cast(long)n[$-1]*N+n[$-2];
      long l=cast(long)m[$-1]*N+m[$-2];
      if(k<l) { r=n; return 0; }
      p = cast(int)(k/l);
    } else {
      p = div32(n[$-3..$],m[$-2..$]);
    }
    assert(p>0 && p<=N);

    if(p==N) {
      p=N-1;
      assert(abscomp(n,absmul(m,[N-1]))>=0);
    }

    Bignum bb=abssub(n,absmul(m,[p]));
    if(bb.sign>=0) {
      r=bb.abs;
    } else {
      p-=1;
      bb=abssub(m,bb.abs);
      assert(bb.sign==1);
      r=bb.abs;
    }

    return p;
  }

  if(m.length==0) throw new ZeroDivException();
  if(n.length<m.length) { rem=n; return []; }

  n=n.dup;
  size_t ln=n.length;
  int[] c;
  c.length=n.length;

  while(true) {
    size_t i=ln-m.length+1;
    assert(1<=i && i<=ln);

    if(abscomp(n[i-1..ln],m)<0) i--;
    if(i==0) {
      rem=n[0..ln];
      break;
    }

    int q; int[] r;
    q=divN(n[i-1..ln],m,r);

    assert(q>0);
    c[i-1]=q;
    if(i==1) { rem=r; break; }

    assert(i>=2);
    if(r) {
      ln=i-1+r.length;
      n[i-1..ln]=r.dup; // n[i-1..ln]とrが重なってるとコピー不可（その場合必ずr.ptr==&n[i-1]でコピー不要な気もするが…）
    } else {
      while(i>1 && n[i-2]==0) i--;
      ln=i-1;
    }

    if(ln<m.length) {
      rem = n[0..ln];
      break;
    }
  }

  return delzero(c);
}

/// 最大公約数（互除法）。ただしgcd(0,0)=0
int[] absgcd(int[] x, int[] y) {
  assert(abschk(x));
  assert(abschk(y));

  if(x.length==0) return y;
  if(y.length==0) return x;

  int c=abscomp(x,y);
  if(c==0) return x;
  if(c<0) {
    int[] t=x; x=y; y=t;
  }

  while(y.length!=0) {
    int[] r=absmod(x,y);
    x=y; y=r;
  }
  return x;
}

/// Bignumの最大公約数（符号に関係なく正）
Bignum biggcd(Bignum x, Bignum y) {
  int[] abs = absgcd(x.abs, y.abs);
  return new Bignum(abs.length ? 1 : 0, abs);
}

/// 絶対値のベキ乗。0**0=1とする。
int[] abspow(int[] n, int[] r) {
  assert(abschk(n) && abschk(r));
  if(n.length==0) {
    return r.length ? [] : [1];
  }

  bool[] a;
  while(r.length>0) {
    int[] u;
    r=absdivmod(r,[2],u);
    assert(u.length==0 || u.length==1 && u[0]==1);
    a ~= u.length!=0;
  }
  int[] ret=[1];
  foreach_reverse(i; a) {
    ret=absmul(ret,ret);
    if(i) ret=absmul(ret,n);
  }
  return ret;
}

/// Bignumのベキ乗(2**-1などは0)
/// Throw: nが0でrが負のときZeroDivException
Bignum bigpow(Bignum n, Bignum r) {
  if(r.sign>=0) {
    int sgn = (n.sign<0 && r.sign>0 && r.abs[0]%2==1) ? -1 : 1; // 結果の符号、nが負でrが奇数なら負
    int[] abs=abspow(n.abs, r.abs);
    return new Bignum(abs.length ? sgn : 0, abs);
  } else {
    if(n.sign==0) {
      throw new ZeroDivException();
    } else if(n.abs.length==1 && n.abs[0]==1) {
      return (n.sign<0 && r.abs[0]%2==1) ? tobig(-1) : tobig(1);
    } else {
      return tobig(0);
    }
  }
}
